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Quasiparticles are an important decoherence mechanism in superconducting qubits, and can be 
described with a complex admittance that is a generalization of the Mattis-Bardeen theory. By 
injecting non-equilibrium quasiparticles with a tunnel junction, we verify qualitatively the expected 
change of the decay rate and frequency in a phase qubit. With their relative change in agreement to 
within 4 % of prediction, the theory can be reliably used to infer quasiparticle density. We describe 
how settling of the decay rate may allow determination of whether qubit energy relaxation is limited 
by non-equilibrium quasiparticles. 

PACS numbers: 74.81.Fa, 03.65.Yz, 74.25. N-, 74.50.+r 
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Superconducting resonators and Josephson qubits are 
remarkable systems for quantum information process- 
ing |J^], with recent experiments demonstrating 3-qubit 
entanglement [2, 13 1 Bell and Leggett-Garg inequalities 
[3, 01 J resonance fluorescence fil, photon jumps [tJ and 
transducers Q, and NOON states j^. For maximum 
coherence, devices are operated at temperatures T well 
below the superconducting transition Tc to suppress ther- 
mal excitations and quasiparticles. Nevertheless, exper- 
iments using resonators and qubits can sense the pres- 
ence of quasiparticles, which arise from non-equilibrium 
sources such as stray infrared radiation, cosmic rays, or 
energy relaxation in materials [13, It is vitally impor- 
tant to completely understand the theory of quasiparticle 
dissipation, determine whether it limits qubit coherence, 
and identify generation mechanisms. 

In this Letter, we describe the effect of non-equilibrium 
quasiparticles with linear response theory [12.] including 
dissipative and dispersive terms, and then experimentally 
test that theory by measuring the response of a supercon- 
ducting phase qubit to the injection of non-equilibrium 
quasiparticles. The decay rate of the qubit is shown qual- 
itatively to depend on the injection and recombination 
rate of quasiparticles. Because the density of quasipar- 
ticles Tiqp is not known from injection, we compare the 
decay rate and frequency shift of the qubit with chang- 
ing riqp to test quantitatively the correct dependence with 
theory. We also show that the characteristic settling time 
to steady state is consistent with theory, and may allow 
an additional measure of quasiparticle density. 

In superconductivity, non-equilibrium quasiparticles at 
energy E can be described using an occupation prob- 
ability f{E). Quasiparticle effects can be expressed in 
a two-fluid manner by considering their total density 
nqp = 2D{Ef) p{E)f{E)dE, where D{Ef)/2 is the 
single-spin density of states, A is the gap, and p{E) = 
E/y/E^ — is the normalized density of quasiparticle 
states. For low dissipation appropriate for qubits, we 



need only consider small occupations riqp ^ ricp, where 
we define n^p = D{Ef)A as the density of Cooper pairs. 
Using ECS theory, non-equilibrium quasiparticles lower 
the superconducting gap as A = Ao(l — nqp/ricp), where 
Ao is the gap without quasiparticles 

We first describe the effect of non-equilibrium quasi- 
particles on Josephson tunneling. For arbitrary occupa- 
tion f{E), the Josephson inductance Lj and admittance 
Yj from Cooper pair tunneling has frequency dependence 
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where (f) is the junction phase difference, $o = h/2e is the 
flux quantum, and Iq is the critical current jl3| . Non- 
equilibrium quasiparticles effectively lower the critical 
current by a factor 1 — 2/ (A) by blocking pair channels. 
Note this reduction is not a function of rtqp , but depends 
on the occupation probability at the gap /(A). For an 
exact calculation, the occupation would be for Andreev 
bound states in the junction at energies slightly below 
the gap |14| : for small tunneling probability, we assume 
this equals the occupation of free quasiparticles at the 
gap /(A). Thermal quasiparticles have an occupation 
fxiE) = 1/[1 + exp(£;/fcr)], which yields the expected 
temperature dependence tanh(A/2fcr). 
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Defining a through the relation /(A) = an, 
find a ~ •^/A/pT'^T^ for thermal quasiparticles at low 
temperature [l5|. For the case of non-equilibrium quasi- 
particles generated by an external pair-breaking source, 
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we numerically compute a ~ 0.12 (nqp/ricp)" 
a ~ 1.2 for typical parameters [l3|. 

The total junction admittance Yj is also affected by 
quasiparticle tunneling [l^ , and is given by 
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where the last term, from Cooper pair tunneling, has the 
critical current Iq = 7rA/2ei?„ re-expressed using the 
gap and junction resistance i?„. The first term gives a 
different phase dependence f -l-cos 0, and has contribution 
from the tunneling of free quasiparticles [riq] and bound 
Andreev states [/(A)]. 

For a junction phase difference ip — the factors from 
/(A) cancel out, and the admittance is 
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This result is equivalent to the normalized conductiv- 
ity a{ijj)/an given by the Mattis-Bardeen theory [16] for 
thermal quasiparticles in the limit kT <C fiw <C A. 
This shows that resistance in a superconductor can be 
modeled as a series of Josephson junctions. 

Quasiparticle damping increases the energy relaxation 
rate in Josephson qubits. For the phase qubit, where 
matrix elements are well approximated by harmonic os- 
cillator values, the relaxation rate is given by (l8| 
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where C is the junction capacitance and Eiq is the qubit 
energy. This result is identical to that found from envi- 
ronmental P{E) theory in the low impedance limit 
and now includes a cos term coming from the interfer- 
ence of electron- and hole-like tunneling events [l9|. 

Quasiparticles also change the imaginary part of Yj , 
which then shifts the qubit frequency by 5EuJh ~ 
—ln\{Yj{Eio/h)}/2C using perturbation theory [18]. For 
the riqp/ncp term in quasiparticle tunneling, the real and 
imaginary parts of Yj are equal, which yields a change 
in angular frequency — Fi/2. Quasiparticles reduce the 
admittance from Cooper pair tunneling by the factor 
1 — (1 + 2a)nqp/ricp due to a change in A and the pair- 
blocking terms. Including the final /(A) term from quasi- 
particle tunneling, the frequency shift is 
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where h = -y/ A/2£'io/7r ~ 0.6 for typical parameters. 
Here, the shift is normalized to the damping since both 
are proportional to nqp. 

The above calculation assumes constant i^, which is 
not valid for a phase qubit where we impose the con- 
straint of constant current bias /. Assuming the standard 
Josephson current-phase relationship, the critical current 
changes both from the gap and quasiparticle excitations 
in the junction, giving SIq/Iq = —(1 + 2a)nqp/ncp. Be- 
cause the qubit frequency scales as Eio oc (/q — /)^^^, a 
change in qubit critical current changes the bias Iq — I 
and results in an additional shift in frequency given by 
5Eio/Eio = {l/4.)SIo/{Io-I)- 



Since the equations for both dissipation and frequency 
shift are proportional to riqp, their ratio can be related 
using qubit parameters. Including both frequency shifts, 
we find fil 
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The latter term dominates because the current bias is 
typically set close to the critical current Iq — I <^ Iq. 

As shown in Fig.[T][a), we experimentally test this the- 
ory by using a phase qubit [20| to measure the effect 
of non-equilibrium quasiparticles injected via a separate 
tunnel junction in the SQUID. The effect of quasiparti- 
cles are modeled by a parallel admittance Yj that changes 
the qubit |1) to |0) state transition frequency and decay 
rate. Qubit measurement produces a state-dependent 
change in the l oop flux that is measured with a SQUID 
readout circuit [2(|. The SQUID is also used to generate 
quasiparticles when driven into the voltage state Vsq > 0. 
With the SQUID shunted by a resistor i?sq = SOfl, the 
SQUID current I^q can be adjusted to produce a voltage 
V^q from ~ 0.6 A/e to above 2A/e, greatly changing the 
generation rate of quasiparticles. 

The experimental time sequence is illustrated in 
Fig. [lib). Quasiparticles are initially generated by apply- 
ing SQUID current to produce T4q ^ 2 A/e for a time Tinj: 
quasiparticles then diffuse to the qubit via the ground 




FIG. 1: (a) Schematic of phase qubit, with critical current 
Jo — 2.0 /lA and capacitance C ~ 1.0 pF, shunted by admit- 
tance Yj coming from non-equilibrium quasiparticles. The 
measurement SQUID, shunted by a resistor i?sq = 30 fi, gen- 
erates quasiparticles when switched into the voltage state, (b) 
Time sequence of experiment, showing initial pulse for quasi- 
particle injection at the SQUID, settling time Tset, flux pulses 
applied to the qubit for frequency or decay time measurement, 
and a final flux measurement by the SQUID. 
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FIG. 2: Effect of quasiparticles on the phase qubit. (a) Plot 
of qubit decay rate Fi (circles) and frequency shift SEio/h 
(triangles) as a function of SQUID injection current /inj, for a 
fixed time of injection rinj = 200 /is and settling Tsct = 300 /is. 
Dashed vertical line indicates SQUID voltage Vsq = 2A/e, 
above which quasiparticle generation increases rapidly, (b) 
Plot of decay rate versus settling time (circles) for fixed in- 
jection current /inj = 20 fiA and time rinj = 200 fis. Quasi- 
particles are observed to recombine on a time scale ~ 300 /is. 
Bottom blue line is theory of Eq. ([7]) with electron-phonon 
coupling [2H To = 400 ns. Theories with additional qubit de- 
cay (top red, To = 200 ns) and non-equilibrium quasiparticles 
of Eq. ([SJ (gray, tq = 400 ns) are also shown, each fit with tq 
and a parameter matching Fi at 10 ms. Qubit number is 2, 
with 7o = 2.03 iiA and C = 1.063 pF. 



plane connection. The injection is turned off for a set- 
tling time Tsct, after which we apply a flux waveform 
to perform a measurement of the qubit resonance fre- 
quency or decay rate ^2d\ . We then ramp /gq to measure 
the flux state for qubit readout, but at a low voltage 
V^q ~ 0.6 A/e so that few quasiparticles are generated. 
The experiment is typically repeated ^ 1000 times to 
produce probabilities. 

We plot in Fig.[5Ja) the decay rate and frequency shift 
of the qubit versus injected SQUID current for a fixed 
settling time. When the current I^q produces a SQUID 
voltage Vsq < 2 A/e (left of dashed vertical line), we ob- 
served almost no frequency shift SEio/h or change in 
the decay rate Fi. However, when the voltage exceeds 
the gap and greatly increases the quasiparticle genera- 
tion rate from pairbreaking, we observe a large increase 
in qubit decay rate and decrease in frequency. 
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FIG. 3: Parametric plot of qubit frequency shift SEio/h and 
decay rate Fi for various injection currents. A linear relation 
is observed, consistent with theoretical predictions that both 
should scale as Wqp. Gray line indicates fit to the data. 

The recombination of quasiparticles is tested by keep- 
ing the injection current and time constant, while varying 
the settling time, as shown in Fig.[5I^b). For short times 
corresponding to the highest quasiparticle densities, we 
see more rapid decay of the qubit. For settling times 
greater than ^ 500 fis the decay rate approaches that 
without injection of quasiparticles. 

Although these data show qualitative agreement with 
expectations, direct quantitative analysis is difficult be- 
cause of uncertainties in predicting Uqp from modeling 
the generation, recombination, and diffusion of the quasi- 
particles. However, an accurate quantitative test of the 
theory can be obtained by comparing the change in decay 
rate with the change in qubit frequency, which are both 
expected to scale as Uqp. 

This comparison is shown in Fig. [3] for the range of 
injection currents plotted in Fig. [2] With both quanti- 
ties proportional to riqp , a linear relation is expected and 
observed in the data. The slope SEio/hTi is thus a quan- 
titative test of the theoretical prediction given in Eq. ([6]) . 
We extracted the slope from this type of plot for 3 de- 
vices and 5 experiments, as summarized in Table I. We 
find that the experimentally measured slope is on aver- 
age only slightly larger (3.6%) than predicted by theory 
and well within experimental uncertainty (14%). 

We have also compared the dissipation and frequency 
shift in superconducting coplanar resonators made from 
aluminum, where quasiparticles are generated by rais- 
ing the temperature. We find the shift in resonance fre- 
quency and inverse quality factor also scale together, with 
a slope that is 0.77 times the predicted valuellSj. The 
difference is probably due to two-level states [23j. 

With confidence that we can accurately extract n^p 
from our experiment, we now analyze the time depen- 
dence of the data plotted in Fig.[2l^b). The decay of 
quasiparticles can readily be calculated for the case of 
small density when they mostly have energy near the gap 
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TABLE I; Qubit parameters showing average cos(j), bias pa- 
rameter (/q — /)//o, operating frequency, slope extracted from 
experiment data as plotted in Fig.O slope predicted by theory 
of Eq. (O, and their ratio. We observe experimental agree- 
ment within experimental uncertainty, as denoted in paren- 
thesis. Qubits 3, 4, and 5 are the same device tested at dif- 
ferent times. Calculation of cos <j) and (Jo — I)/Io accounts for 
the inductive shunt in the flux-biased phase qubit. 
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(GHz) [4§f]Exp 



Exp 
Thy 



1 0.12(3) 0.042(5) 6.523 

2 0.17(4) 0.066(6) 6.743 

3 0.19(4) 0.071(3) 6.413 

4 0.33(5) 0.130(3) 7.375 

5 0.19(4) 0.071(4) 6.413 



-1.42(9) -1.34(19) 1.05(16) 

-1.30(3) -1.03(13) 1.22(15) 

-1.13(6) -1.03(09) 1.09(12) 

-0.67(2) -0.80(14) 0.83(14) 

-1.04(6) -1.04(10) 0.99(11) 



[l^ . In this limit, the integral of Eq. (7) in Ref. gives 
a recombination rate — (21.8/To)(nqp/ncp), where 
To ~ 400 ns is the characteristic electron-phonon coupling 
for aluminum |2l| . If rqp is the injection rate, the differ- 
ential equation dnqp/dt — — 2r''nqp 
to give 



rqp can be solved 
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with equilibrium values (r'')'^i = [(43.6/ro)(rqp/ncp)] 
and (nqp)°i/ncp = (To/43.6)(r'')''i , where to is the in- 
tegration constant. After quasiparticle injection, Eq. ([8]) 
describes a decay that initially follows Eq. (O, but levels 
off to a steady state value (nqp)*"^ after time ~ l/{T^y^. 

We plot in Fig.[2)Jb) the prediction for the case of no 
equilibrium quasiparticles Tqp = [Eqs. ^ and ([7])], 
which depends only on tq. To account for the constant 
rate at times > 1 ms, we consider two hypotheses. One 
is the qubit has an additional dissipation mechanism, so 
that the net decay rate is the sum of the above predic- 
tion with a constant offset. The second assumes constant 
injection of non-equilibrium quasiparticles, as given by 
Eq. ([8]). We plot these two additional predictions, which 
show somewhat different dependencies with time. 

For settling times < 300 /j,s, we believe the departure 
from theory is due to diffusion effects that are not in- 
cluded in this simple model Numerical simulations 
[l^ show such deviations are reasonable given that the 
injection and qubit junctions are well separated. Note the 
small dip in Fi at time ^ 1 ms is similar in magnitude to 
prior temperature measurements (TT| . Although both hy- 
potheses can roughly explain the data for times > 300 /is, 
the scenario that fits a bit better is when non-equilibrium 
quasiparticles are the dominant decay mechanism. 

Fluctuations in the quasiparticle density cause the 
qubit frequency to jitter, producing dephasing [22|. Us- 
ing a Ramsey fringe protocol, measurements were made 
of the decay time T2 at different quasiparticle injection 



currents, where we saw a decrease in T2 at high quasi- 
particle densities. After correcting for the decrease in T2 
from a change in Ti using 1/12 = 1 /2Ti + 1 /T^, we found 
that T(j, was unchanged to within experimental error. We 
conclude that other dephasing mechanisms are dominant 
in the phase qubit. 

In conclusion, we have used the theory of quasiparti- 
cle admittance to predict the qubit frequency shift and 
energy decay from non-equilibrium quasiparticles. When 
injecting quasiparticles from a nearby junction, we find 
good qualitative agreement with theory. Good quantita- 
tive agreement was observed between the relative change 
of qubit decay and frequency with changing quasiparticle 
density. Monitoring the time dependence of qubit decay 
after injection should allow future experiments to pos- 
itively identify if non-equilibrium quasiparticles are the 
limiting decay mechanism in superconducting qubits. 
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Supplementary material is presented for the paper "Energy decay and frequency shift of a super- 
conducting qubit from non-equilibrium quasiparticles". First, we discuss quasiparticle data for a 
superconducting coplanar resonator. We then document how the energy dependence of the occupa- 
tion can be calculated numerically. We calculate analytically the total quasiparticle decay rate and 
time dependence with both a bulk model and, numerically, one including diffusion. We also com- 
pute the quasiparticle dependence of the gap, occupation probability, current-phase relationship, 
and how the frequency shift and dissipation are related. Finally, we calculate the Josephson current 
for non-equilibrium quasiparticles. 
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PACS numbers: 74.81. Fa, 03.65.Yz, 74.25.Nf, 74.50.+r 

We are interested in how non-equilibrium quasiparti- 
cles affect the properties of a Josephson junction in qubit 
devices. Since we are concerned with very low temper- 
ature operation, we consider phonon temperatures suf- 
ficiently low that no quasiparticles exist due to thermal 
generation. The non-equilibrium quasiparticles are gen- 
erated with an unknown mechanism, and then relax their 
energy via the emission of phonons. The quasiparticles 
typically have energy E close to the gap, from which 
they eventually decay via recombination. From electron- 
phonon physics, we know that the qusiaparticles relax to 
energies very near the gap, as calculated in Ref. 

The factor of 2 in the definition of riqp comes from in- 
tegration only over positive energy, whereas excitations 
arise from both electron states above and below the Fermi 
energy. Note this integral already contains the two pos- 
sible spin states in the definition of D{Ef). 



QUASIPARTICLES IN COPLANAR 
RESONATORS 



The relationship between quasiparticle damping and 
frequency shift may also be tested in superconducting 
resonators. Here microwave transmission is measured to 
extract the resonance frequency / and the quality fac- 
tor Q, with the quasiparticle density changed by simply 
increasing the temperature [ij. As shown in Fig.[Tl we 
find that with increasing quasiparticle density, dissipa- 
tion increases and the resonance frequency decreases, in 
a similar manner as for junctions. 

To compare with theory, we use solutions to the 
Mattis-Bardeen conductivity that is valid for the regime 
kT ^ hu! <ti A 0, 0] ■ These results can be expressed in 
terms of an admittance function discussed in the main 
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FIG. 1: Parametric plot of fractional frequency shift Sf/f 
versus dissipation 1/Q, for various quasiparticle occupations 
riqp changed by varying the the sample temperature. The de- 
vice is a aluminum coplanar resonator fabricated on sapphire. 
The slope of the data is in good agreement with predictions 
(gray line) at temperatures above 250 mK. Deviations at low 
temperature are believed to come from the two-level states. 



article, but with a modification of the the l-|-i term 
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dn — 2\/2x/tt sa\\\{x)KQ{x) , 
di — \J2-Kx exp(— a;)/o(a;) , 



(1) 
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where x = hjjj /2kT . The theoretical prediction, indicated 
by the gray line in Fig.[Tl is in reasonable agreement with 
the data at high temperatures. The experimental slope 
at temperatures above 200 mK is 0.77 times that given 
by theory. The deviation is largest below about 120 mK, 
and believed to arise from two-level states that are not 
included in this model. 
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NUMERICAL SOLUTION OF QUASIPARTICLE 
RECOMBINATION 



For numerical computations of non-equilibrium quasi- 
particle density from relaxation and recombination 0], 
the integrals over energy have to be put into discrete 
form. For a binning size given by de in energy, we define 
the number of excitations in bin i as 

Ui = dep{ei)f{ei) . (4) 

Using this definition, the total quasiparticle density is 



(5) 



The scattering and recombination rates of Eqs. (6) and 
(7) of Ref. can then be expressed in a discrete form 
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where the phonon occupation factor is Np{E) = 
1 / \ exp{— E / kTp) — 1| and the bracketed terms arc the 
G factors. We have also assumed small occupation, so 
that in Eq. © we use 1 — / ^ 1. 

The coupled differential equations for the change in the 
excitation number are 



dt 



J2Gt,n,~J2^ + S,,)G:^n,n, , (10) 



where 5ij is the Kronecker delta and accounts for the 
annihilation of 2 quasiparticles when in the same bin (i = 
i). With the physics expressed in matrix form, a solution 
can be readily solved numerically. 



QUASIPARTICLE DECAY 

The physics of quasiparticle relaxation and recom- 
bination was discussed in Ref. 4]. Although the arti- 
cle described solutions for the non-equilibrium occupa- 
tion f{E) using numerical methods, quasiparticle de- 
cay physics can be understood in the case of low den- 
sity where they are mostly occupied at the gap. The 



electron-electron recombination rate of a single quasipar- 
ticle, starting from Eq. (7) of Ref. can be well approx- 
imated using 
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where we have used the BCS result ^jkT^ — 1.76. Here, 
D[Ef)/2 is the single-spin density of states, and we define 
the Cooper pair density Ticp = D{Ep)A. 

The time dependence of the quasiparticle density can 
be understood via the rate equation 
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where a recombination event removes 2 quasiparticles, 
is the single particle quasiparticle injection rate. The 
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second equation is for the normalized quasiparticle den- 
sity, and has a recombination rate that is proportional 
to Tiqp because of the two-body electron-electron interac- 
tion. 

The equilibrium quasiparticle density is given by set- 
ting driqp/dt = 0, yielding density and recombination 
rates 
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The first equation is close to what was found numerically 
in Ref. Q . The second is given by the geometric mean of 
the normalized injection and the characteristic electron- 
electron interaction rates. 

We compared the results of this simple calculation with 
numerical solutions for a range of injection rates and 
found excellent agreement for riqp/ricp ^ 0.001. Even 
at large density riqp/ricp = 0.1, Eq. (IT7|) is a reasonable 
approximation as its prediction is only 40 % larger than 
that obtained via numerics. 

For no injection of quasiparticles rqp = 0, the differen- 
tial equation can be integrated to give 
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where t is the time and to is an integration constant, 
which is approximately the time at which the quasiparti- 
cles start to cool. The solution to the differential equation 
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for a finite injection rate is 

nqp = {n^pr coth[ {T^nt -to)], (20) 

where the coth term is replaced by tanh if the quasiparti- 
cle density increases with time. At short times the term 
cothr'^t = which then gives Eq. (|19|) and a time 

dependence that scales only with the electron-phonon 
coupling time tq. There is a relatively sharp crossover 
to the long time behavior where the quasiparticle den- 
sity riq^ is constant with time. The crossover time is 
given by l/^Vfi. 

The inverse of the crossover time thus gives the equi- 
librium recombination rate (r'')'^i, which is related to the 
density using Eq. (1171) and the parameter tq. Comparing 
this density with that found from the qubit decay rate 
allows one to determine whether quasiparticles are the 
limiting decay mechanism for the qubit. 



QUASIPARTICLE DECAY WITH DIFFUSION 

The analysis in the last section assumes a bulk (uni- 
form) model where there is no diffusion of quasiparticles. 
Here, we describe a numerical solution for quasiparticle 
decay including relaxation, recombination, and diffusion 
using the simple geometry of a thin superconducting disk 
of radius 5 mm. We use constant quasiparticle injection 
throughout the disk and a large injection pulse into the 
center of the disk at time t = —200 /is to t = 0. Because 
diffusion depends on the quasiparticle energy, the calcu- 
lation keeps track of the occupation probability for both 
the radius and energy variables. 

In Fig. [2] we plot quasiparticle density versus settling 
time in a manner similar to that in the main paper, but 
for 4 radii. We find differing behavior depending on the 
ratio of the radius with the diffusion length ~ 1 mm, as 
computed for e-e diffusion in Fig. 3 of Ref . 4] . For small 
radii, we see a dependence on time that matches closely 
with the bulk theory, as described in the main paper. For 
a radius much larger than the diffusion length, the quasi- 
particle density does not change. For the radius close to 
the diffusion length, we observed behavior between the 
two limits - a reduced peak density but a relaxation to 
the steady state value that has a similar time scale than 
for a small radius. 

We note that the actual qubit device has interruptions 
in the ground plane due to the device geometry, so that 
this computation will not exactly match the experimental 
data. However, the model mimics the time dependence 
of the data fairly well above 10 /is, so it is reasonable to 
compare to the simple bulk analysis for the behavior at 
long times > 250 ^s. 
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FIG. 2: Plot of quasiparticle density versus settling time for a 
superconducting disk of radius 5 mm with quasiparticle injec- 
tion at the center. Points are for the simulation at four radii, r 
= (black), 0.5mm (red), 1.25 mm (green), and 5mm (blue). 
Large changes are observed for small radii, which show a time 
dependence close to that predicted by the full theory of Eq. (9) 
of main article (black line) and for the zero background of Eq. 
(8) (cyan line). At large radii much greater than the charac- 
teristic diffusion length of ~ 1 mm, no change in quasiparticle 
density is seen. The simulation for radius 1.25 mm is in rea- 
sonable agreement with experimental data. 



DEPENDENCE OF GAP ON QUASIPARTICLES 

The change in the superconducting gap A with quasi- 
particles can be calculated starting from the BCS gap 
equation, but assuming a small non-equilibrium popula- 
tion fiE) 

A =DiEf )v£'' dEp^il - 2/) , (21) 

where V is the attraction potential and 0£) is the Debye 
energy. Solving for the gap, one finds 

~A„(l _ ^3E) , (26) 

where Aq is the normal expression for the BCS gap with 
no quasiparticles. 
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FIG. 3: Plot of a = /(A)/(nqp/ncp) versus riqp/ncp, found 
from numerical simulation. Based on the values of riqp/ncp 
in Fig. 2b in the paper, we find a ~ 1.2. For a wide range of 
injection rates, a can be well approximated by the power law 
a~0.12(nqp/ncp)-°■"^ 



(riqp /ricp) /S./huj. The Andreev bound states of the 
quasiparticles have an inductance, so the AC Josephson 
relation can then be used to find the current 

/abs(0) - -"^ C l^WYABsm d<P (29) 
= /o[<^ + Sin A) . (30) 

Note that this term increases the critical current, and if 
one replaces — ^ sin it cancels out the decrease coming 
from Josephson tunneling. Since many experiments have 
shown the temperature dependence of the current-phase 
relation is given by only Josephson tunneling, we do not 
use this Andreev bound state term in our calculations. 
In addition, our data is not consistent with including this 
term since it has the effect of reducing 2a in Eq. (|28l) to 
a value below a. 



NUMERICAL DETERMINATION OF 
OCCUPATION PARAMETER 

To determine the effect of non-equilibrium quasiparti- 
cles, both the quasiparticle density riqp and the occupa- 
tion probability at the gap /(A) must be calculated. We 
plot in Fig. [3] the quantity a — /(A)/(nqp/ncp) versus 
nqp/ricp obtained from numerical computations for wide 
range of injection rates. We find the results are well ap- 
proximated by a line on the log-log plot, implying that 
the dependence can be well approximated by the power- 
law formula a ~ 0.12 {nqp/nc-pY 



-0.173 



CURRENT-PHASE RELATIONSHIP WITH 
QUASIPARTICLES 

The current-phase relationship from Josephson tunnel- 
ing is given by 



I{4>) = Iq sin ( 
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[1 - 2/(A)] , (27) 



where Iq = tt Aq /2e_R„ and the dependence of A on quasi- 
particles is now explicitly shown. To first order, the frac- 
tional change in critical current is 
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An interesting question is whether the quasiparticle 
tunneling terms should also be included in the current- 
phase relation. For the junction current to only be a func- 
tion of phase, it must arise for a purely inductive com- 
ponent of the junction admittance, which corresponds 
to terms with a frequency dependence that scales as 
1/iuj. Tunneling of free quasiparticles should not be in- 
cluded since it has an additional frequency dependence 



RELATING DISSIPATION AND THE 
FREQUENCY CHANGE 

Since both dissipation and the fractional critical- 
current change are proportional to riqp, the magnitude 
of these effects are related. The fractional change in the 
qubit resonance frequency Eio/h can be calculated know- 
ing that its dominant scaling is Eio/h cx (/g — lY^'^, 
where / is the qubit bias current, giving 
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Quasiparticle dissipation can be likewise written in 
terms of the quasiparticle density, provided we first re- 
express the capacitance C into qubit parameters. Using 
the Josephson inductance Ljq = ^o/2t:Iq, the qubit res- 
onance frequency is given by 
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We thus calculate the decay rate of the qubit 
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By taking the ratio of Eqs. ((33|) and ((38t . the quasipar- 
ticle densities cancel out, and we can relate the dissipa- 
tion to the frequency shift 
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JOSEPHSON EFFECT FOR ARBITRARY 
QUASIPARTICLE OCCUPATION 

Here we calculate the effect of a non-equilibrium pop- 
ulation of quasiparticle states on Josephson tunneling, as 
appropriate for qubit devices. The current proportional 
to cos 5 is also evaluated for the case of gaps that are not 
equal. The results are readily obtained using standard 
second-order perturbation theory and simple integration 
of intermediate formulas. 

The work expands on Ref. 0, which calculated the 
Josephson effect at zero temperature. In the paper, the 
section on Josephson tunneling is the starting point of 
this calculation. 

The Josephson effect is derived by calculating the 
second-order change in energy to a superconducting state 
from a tunnel junction. The tunneling Hamiltonian in 
second-order perturbation theory is given by 



H. 



(2) 



Ht — Ht 
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where is the energy of the intermediate state i. Be- 
cause the terms in Ht have both and 7 operators, 
the second-order Hamiltonian has terms that transfers 
charge across the junction but does not change the su- 
perconducting state, thus giving a change in the energy of 
the state. This differs from first-order tunneling theory, 
which produces current only through the real creation of 
quasiparticles. 

Because Ht has terms that transfer charge in both 
directions, HtHt will produce terms which transfer two 



electrons to the right, two to the left, and with no net 
transfer. With no transfer, a calculation of the second- 
order energy gives a constant value, which has no physical 
effect. We first calculate terms for the transfer of two 
electrons to the right from (^t+ + ^t-)(Ht+ + Ht~)- 
Nonzero expectation values are obtained for only two out 
of the four terms, as given by 
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12 (cLC_i)(4cLfl,) + {c^LCL)i(^-RCR) 



(43) 

The pairs of electron creation and annihilation opera- 
tors can be computed, giving 



CfcC_fc = (m7o + we"^7j)(u7i - ve"''jl) 

^- W^(-7o7^ -f 7|7i) , 
C-fcCfe -> wwe*"*(-7($7o + -/^-/l) , 
4cLfc = (w7(l + we""^7i)(u7j - ve~"^-/„) 
^ uwe-*"*(-7oSo +7i7j) , 



(44) 
(45) 
(46) 
(47) 
(48) 
(49) 



where we have only included pairs of quasiparticle opera- 
tors 7 that leaves the superconducting state unchanged, 
as needed for a calculation of the energy change from 
tunneling. Inserting these operators into Eq. (j43p and 
defining the phase difference 5 = 0l — we find 



(-'7lo71o + 7li7Li)(-7™7fl,o + ^riIri) + {^llolha + 1 Lilli){-1 Ralm + IriIri) 



(50) 
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El — En 
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where we have computed the intermediate energy ti using a positive (negative) energy E for the creation (annihilation) 
of a quasiparticle. The quantity uv = IS./2E describes the amplitude for the virtual quasiparticle to be both electron- 
and hole-like, which allows a net transfer of charge by two electrons. 
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We now change the sum to an integral over electron states according to 

/OO /"OO /'OO /'OO 

d^L Nor / dU = 2NoL / Pl dE^ 2Nor / pR cLEr , (52) 



where A^'o is the normal density of states, and p = E/y/E^ — is the (normalized) superconducting density of states. 

By describing the superconducting state with an occupation probability of quasiparticles / = f{E), the quasiparticle 
operators for the creation then destruction of a quasiparticle is weighted by 1 — /, while the process of destruction 
then creation is weighted by /. The tunneling Hamiltonian is then given by 

h!^> = \t\^e''NoLNoR 2 / pLdEL 2 / pRdER 2G , (53) 

JAl JAr ^El ZEr 

(l-/L)(l~/fl) hfR a-fL)fR fdl-fn) , . 

El+Er ~El-Er EL-ER + ie -El + Er + Ie ^ ' 



^ — Jl - Jr + JlIr , /i/fl. , Ir - JlIr. , Jl - /l/r 



El+Er El+Er El-Er + ie -El + Er 

^ - Jl - Jr , „ ^ fL 



(55) 



p ^ p . ^ P P + ^Ah + fR- 2fLfR)5{EL - Er) (56) 

i^L + J^R i^L — ^R 

1^ + P ^^^^rp^'^^ + + fR- 2fLfR.)S{EL - Er) , (57) 

J^L + t^R t^L ^ ^R 



where G comes from quasiparticle operators (bracket 
terms in Eq. (|5ip ) after removing a factor of 2 because 
of the pair of states and 1. We take e — >■ 0+, and the 
integration over the zero of energy in the denominator is 
performed using l/{x + it) — P(l/x) + in5{x), where P 
is the principle part and 5{x) is the Dirac (5-function. 

The total second-order Hamiltonian for the tunneling 
of 2 electrons in both directions is 

= ll^ + h.c. (59) 
= 2Re{i4^} , (60) 



where h.c. is the Hermitian conjugate. 

This result can be expressed in more physical terms by 
noting that the junction resistance can be written as 

1 47rp2 

= \t\'NoRNoL ■ (61) 

Rn h 

In addition, The Josephson tunneling current is given by 

' - h dS ■ ^ ' 

Combining all of these equations, the total Josephson 
current is given by the integrals 
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Wc note that the sin 5 term in Eq. (l63l) corresponds to tegration over El 
Eq. (22) of the Ambegaokar-BaratofF calculation Q. 



We evaluate these integrals by first considering, with- roo ^ ^ 

out loss of ge nerality, that A^ < Ar. Using pA/E = = '^^ / 2 _ ^^^^ 
A/VE^ - A2 and y = Er/Al > 1, we compute that the arccosh ^ 

temperature independent term 1/(El + Er) gives for in- = ^ — (66) 
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The remaining integration over En gives 0] 

An ALarccosh(£'i^/AL) 



IlLR 



dEii 



TT — EUipticK 
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(for Al = Ar = A), 



(67) 
(68) 

(69) 



be evaluated 

hRL = -ttMAl) ALEllipticK[l - [AL/ARf] (80) 



-A 2/i(A) 



(81) 



where in the last equation we have taken the limit A^ — > 
Ar = A. 

For the case of equal gaps A, the first two integrals 
give a Josephson current with a sin S dependence 



where the last equation uses EllipticK(O) — Tr/2. From 
numerical integration, we have found that Eq. (j68p is only 
approximate for Al ^ Ar. 

For the next term in Eq. (|63p that has the principle 
part of El/{E\ — Ej^), we first integrate over Ei^. With 
the assumption A^ < A^, the integral always passes 
across the pole at Er giving 
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We thus find no contribution for /r in the total integral. 

We next compute the Er integral for the Er/{E\ — 
Ej^) term. We define w = El/Ar, and note that the 
integration over Er depends on whether El is greater or 
less than Ar 
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This result implies that when integrating over E^, no 
contribution comes from Ar to infinity, so the full inte- 
gral is 
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In the limit where A^ is close to Ar such that is 
constant over the region of integration, the integral can 
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2 eRn 

The last equation assumes a thermal population of 
quasiparticles given by f{E) = 1/[1 -I- cxp(i?/fcr)], 
which yields the Ambegaokar-Baratoff formula [6] for the 
Josephson current. 

The Josephson current can also be calculated for an 
arbitrary quasiparticle occupation under the assumption 
that the difference of the gaps Ar — Al is much larger 
than the typical width of the quasiparticle distribution. 
The contribution from the principle part for the /j^, term 
is zero, as discussed before. The contribution from /l 
(the lower gap side of the junction) is given by Eq. (|79p . 
Noting that Jl is peaked at = A^,, we find the current 
change from quasiparticles is given by 

cAb a 



2 sin S 
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sinS ArAl riqpL 
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h dEL 



(85) 
(86) 



Note that these equations have a contribution from quasi- 
particle occupation only from the left electrode, which 
has lower gap. This makes sense since a more exact the- 
ory of Andreev bound states has suppression of the crit- 
ical current from occupied states in the gap, which has 
to have energy below that of the lowest gap. 

For the cos S term, we first note that the current di- 
verges logarithmically for Al — >■ Ar. Assuming the 
quaisparticle density is constant with / = .fi + Jr — JlIr, 
numerical integration gives the approximate formula 

For the case where the gaps greatly differ, and in the 
limit discussed in the previous paragraph, the current is 



2cos(S .AlAr 
—Pl{Ar) I PRdEfR 
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nqpR 



eRn yA2 



Al NorA 



R 



(89) 
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which has a form and magnitude similar to the thermal 
current. 

Note that for the cos J current, quasiparticles con- 
tribute from the higher gap side of the junction. This 
contrasts the behavior of the sin 5 current, which has con- 
tribution from the superconducting electrode with lower 
gap. 
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